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Birefringent crystals are used to control and modify po- 
larized light in a wide variety of applications from the 
millimeter-wave to the ultraviolet. Jones- or Mueller- 
matrix formalisms, which couple incoming and outgo- 
ing polarization states represented as vectors, offer use- 
ful characterizations of such devices. Idealized Jones and 
Mueller matrices for most devices arc found in stan- 
dard textbooks on the subject [1-3]. These ideal matrices 
are adequate in many applications; however, a general 
method for computing Jones and Mueller matrices for 
non-ideal devices including the effects of multiple reflec- 
tions between devices is important where high precision 
and tight control of systematic errors is necessary. 

Real devices significantly diverge from ideal behav- 
ior when used across a wide frequency band or a range 
of angles of incidence. In addition, anti-reflection coat- 
ings that are important in creating high-throughput, 
bircfringcnt-crystal devices are themselves sources of 
non-ideal behavior. And often wave plates are made of 
multiple layers of birefringent materials oriented relative 
to one another in such a way as to increase the frequency 
bandwidth [4]. 

Matrix methods for treating stratified media have a 
long history. Abeles developed a 2 x 2 characteristic ma- 
trix method for fast computation of transmission and re- 
flection for stratified isotropic media [5]. In Abeles' treat- 
ment separate characteristic matrices were required for 
s- and p-polarized plane waves. See also [6, 7] for treat- 
ments of isotropic stratified media. A number of authors 
have expanded to using 4x4 matrices to treat birefrin- 
gent materials that introduce mixing between orthogo- 
nal polarizations. Berreman and others developed matrix 
methods written in terms of first-order differential equa- 
tions [8]. Yeh derived a matrix method coupling propa- 
gating electromagnetic modes [9]. 

We present an alternative derivation of a 4 x 4 ma- 
trix that directly couples total-field components at the 
interfaces between layers in a stratified, possibly bire- 



fringent, medium. The matrix method developed allows 
for full calculation of the polarization state of both trans- 
mitted and reflected waves at any frequency and angle of 
incidence of a stack of any number of isotropic and bire- 
fringent crystals, assumed to vary only in the z direction 
and be infinite in the x- and y-directions. 

This treatment allows direct connection with the Jones 
and Mueller matrices, as shown below. Explicit formulas 
for a case of great practical importance, the matrix of a 
uniaxial crystal with its optic axis in the x-y plane, are 
derived. These results were first published in the author's 
doctoral dissertation [10] with particular application to 
millimeter-wave polarimetry. 

1. Parametrizations of Polarization 

The polarization state of a plane electromagnetic wave 
can be represented in multiple equivalent ways. The 
two-element Jones vector, {\E X \ exp(iS x ), \E y \ exp(iS y )), 
gives the amplitude and phase of two orthogonal, time- 
harmonic electric fields as complex numbers. One can al- 
ternatively parametrize the polarization state of a plane 
wave using the Stokes parameters, defined as 

P = (EEl)=Ia I + Qa Q + U<T U + Va v , (1) 

where the a l are the Pauli matrices 1 , E is the complex 
electric field vector, E^ is its complex conjugate, and an- 
gled brackets represent averaging over a time that is long 
compared with the period of the electromagnetic wave, 
but short compared with time scales of the measurement. 

1 Note that this ordering of the Pauli matrices is different from 
that often used in the literature. 

ai = ( 1 ) CT « = ( -1 ) 
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The four Stokes parameters can be gathered into a vec- 
tor (7, Q, U, V) with 4x4, real- valued Mueller matrices 
relating incoming and outgoing Stokes vectors. Given a 
Jones matrix, the expression 



(3) 



where i and j are in the set {I,Q,U,V}, can be used 
to calculate the corresponding Mueller matrix [11]. It 
should be noted that Mueller matrices can be band- 
averaged element-by-element [12] making the matrix 
method described applicable to experiments with finite 
bandwidth in a straightforward way. 

2. The generalized transfer matrix 

We seek to relate the total, tangential, complex D and 
H field components at the two interfaces of a system 
that is infinite in the x and y directions and varies 
only in the z direction. The plane of incidence is as- 
sumed to be the x-z plane. We will connect the vec- 
tor A = (D,W H,W D,M -HW) on the two interfaces. 
The great advantage of the generalized transfer matrix 
is that it directly relates the field amplitudes at the 
two interfaces, taking into account multiple reflections 
and allowing the optical properties of a system involv- 
ing arbitrarily-many layers to be calculated by simple 
matrix multiplication. 

Once the generalized transfer matrix is computed for a 
system, it is straightforward to calculate the amplitude 
transmission, reflection, and absorption coefficients for 
a wave incident upon the system from one isotropic di- 
electric of index m and transmitted through to another 
isotropic dielectric of index 713 . This situation is depicted 
in Figure 1, where the specific case of a birefringent crys- 
tal is shown; however, we are presently only concerned 
with the fields on the outsides of the crystal. For our 
present purposes, one could replace the anisotropic crys- 
tal with any system of isotropic and birefringent crystals 
represented by a total generalized transfer matrix, T. If 
we denote reflected and transmitted waves by subscripts 
r and t, respectively, and s- and p-polarizations with su- 
perscripts, the fields at the interfaces are related via the 
generalized transfer matrix as 
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The plus and minus signs between incident and re- 
flected electric field components on the left-hand side of 
this equation are dictated by the boundary conditions. 
When the z component of k changes sign for the inci- 
dent versus reflected waves, the z components of E and 
H follow suit according to H = nk x E. The only way 
to simultaneously satisfy the boundary condition in the 
x components of those fields is with a 180° phase shift in 



those components. This introduces a minus sign in the 
first and last elements of A. 

Also, note that while T relates components of D, 
Equation 4 is written in terms of E, which is simply equal 
to D/n 2 in the isotropic media assumed to surround the 
system. This is done to facilitate the calculation of the 
Jones matrix below. 

A. Jones and Mueller matrices 

The above system of equations can be solved to give E£, 
E* , Ef , and Ef in terms of Ef and Ef . Our goal is to write 
these relations in the form of Jones matrices for trans- 
mission and reflection, which can then be translated to 
Mueller matrices using Equation 3. Denoting the com- 
ponents of the 4x4 generalized transfer matrix as t tJ , 
we can simplify our expressions by defining 
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In terms of these constants, the transmitted-wave ampli- 
tudes are 
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The Jones matrices for the reflected and transmitted 
waves can then be transformed into Mueller matrices 
using Equation 3. 

3. Matrix of a uniaxial crystal 

We now derive the transfer matrix for a uniaxial, bire- 
fringent crystal with its optic axis parallel to the faces of 
the crystal in the x-y plane. For full treatments of wave 
propagation in birefringent media and refraction at an 
interface between birefringent media see [6, 13, 14]. 
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Fig. 1. Geometry of the rays used to calculate the generalized transfer matrix for the uniaxial crystal. The two 
polarization states defined by the plane of incidence, the s- and p-waves, are mixed inside the uniaxial crystal into 
the ordinary and extraordinary waves. 



A. Electromagnetic waves in a uniaxial crystal 

Anisotropic, non-magnetic, linear dielectrics can be char- 
acterized via their dielectric tensor, Cki , which transforms 
the electric field vector in the material to the electric dis- 
placement. In the special case considered here of a uniax- 
ial crystal with its optic axis tangent to interfaces I and 
II and rotated at an arbitrary angle, to the x-axis, 

/ n 2 e \ 
e' = R( X )[ n 2 )R (-*) , (9) 
V nl J 

where R{x) is the rotation matrix by \ about the z axis 

(cos x — sin \ \ 
sin X cos X (10) 
1/ 

The above matrix holds for a uniaxial crystal with its op- 
tic axis pointed in the direction of £ = (cosx,sinx,0). 
The crystal is assumed to have ordinary and extraordi- 
nary indices of refraction of n and n e . The inverse to e' 
transforms D to E. 

One consequence of the anisotropy of the crystal is 
that the electric field vector, E, and the electric dis- 
placement, D = sE, of an EM wave no longer point 
in the same direction. Both E and D remain perpen- 
dicular to the magnetic field, H , and the direction of 
energy transport is still given by the Poynting Vector, 
S = E x H; however, the wave propagates in the direc- 
tion of D x H, given by the normal direction to planes 



of constant phase. It is therefore the triplet of vectors 
k, D, and H that are mutually orthogonal. This makes 
it useful to use boundary conditions on the components 
of D at the interface between a birefringent crystal and 
an isotropic dielectric. This is especially true since the 
magnitudes of D and H arc related simply by 

\H\ = -\D\. (11) 

n 

Given k and either £ or D, the other vectors for the 
wave can be calculated through Equations 9 and 10, as 
well as the fact that 

H=-kxD. (12) 
n 

Here n is a real- valued, angle-dependent refractive index 
which gives the speed of a ray as it traverses the medium, 
v p = c/n, that will be calculated below. 

Another consequence of the anisotropy of the mate- 
rial is that a linearly-polarized plane wave at a dielectric 
boundary gets refracted into two separate plane waves 
polarized along preferred directions in the crystal unless 
the incoming plane wave happens to be polarized along 
one of the principle axes of the crystal. These two waves 
travel at different speeds. For a uniaxial crystal, one of 
these waves, the ordinary ray, always travels at the same 
speed, given by 

v' = v = c/n , (13) 
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while the extraordinary ray has a velocity which depends 
upon the angle ip between the wave and the optical axis 
and is given by 



(yl cos 2 ip + v\ sin 2 ip) 

1/2 



1/2 



( cos 2 ip i sin 2 ip \ 
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An EM wave which is incident at angle 9\ upon a sur- 
face of the uniaxial crystal from an isotropic dielectric of 
index m will then refract by two separate angles, which 
are each individually given by the familiar Snell's Law 



of propagation k and the direction of the optical axis 
Thus 



= a k x £ 




(20) 



Similarly, the extraordinary ray must have an electric 
displacement that is perpendicular to its direction of 
propagation and to the direction of vibration of the or- 
dinary ray. It is thus 



n Q sin 9' = n" sin 6" = n\ sin 9i . 



(15) 



The two rays remain in the plane of incidence. The or- 
dinary ray is always at the same angle relative to the 
angle of incidence. The equation for the extraordinary 
ray is complicated by the fact that the index of refrac- 
tion depends on the angle through the material, which 
depends on the angle of incidence and orientation of 
the crystal, x, in a nontrivial manner. Let the faces of 
the uniaxial crystal lie in the x-y plane. The extraordi- 
nary ray that propagates through the material at angle 
6," yet to be determined, has unit propagation vector 
k = sin 9" x + cos 9" z . The angle ip, defined in Equation 
14, between these two is then cos tp = k" •£ = sin 9" cos %• 
Thus Equation 14 becomes 



ip , sin ip 
I ^2 



-1/2 



n a n e [(n 2 - n 2 ) sin 2 9" cos 2 \ + n l] 1/2 ■ 



and Equation 15 becomes 
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Solving this equation for sin 9" yields 
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(18) 



Equations 15 and 18 give the angles of the ordinary 
and extraordinary rays for a given angle of incidence and 
rotation of the HWP. These rays both travel in the x- 
z plane and have associated unit-normal vectors k and 
k , which are 



k = (sin 



6»',0 



l,cos6"); k" = (sin 6>", 0, cos 6>") (19) 



The directions of vibration of the electric displacement 
vector, or in other words the polarizations of the two 
refracted rays, can also be calculated. The ordinary ray, 
traveling through the material at angle 9' , must have its 
electric displacement perpendicular to both the direction 
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— cos x cos 9' cos 9" 
a" \ - sin x [sin 9' sin 9" + cos 9' cos 9"] 
cos x cos sm ®" 

The constants a' and a" normalize these vectors. 

B. Relations between field components 

The geometry under consideration is shown in Figure 1. 
A plane wave of frequency v is incident on a bircfrin- 
gent dielectric. The plane of incidence is the x-z plane. 
We wish to derive relationships between tangential D 
and H at both of the interfaces, labeled I and II. In the 
absence of free charges on the boundary, the tangential 
components of both D and H are continuous. We will 
thus consider the relationship between the vector (D^, 
Ky, D^, -H^) on the two boundaries. In the calculations 
below, subscripts will denote incoming, transmitted, and 
reflected field amplitudes as lower-case i, t, and r, respec- 
tively, along with a roman-numeral I or II denoting the 
interface concerned. Subscripts s and p will denote waves 
with electric field vector perpendicular to the plane of in- 
cidence (s-polarized wave) and in the plane of incidence 
(p-polarized wave) in the surrounding isotropic medium. 
In the uniaxial crystal ' or " will denote the ordinary 
or extraordinary rays. As an example, E"j is the electric 
field amplitude of the extraordinary ray that is transmit- 
ted from interface I. A subscript of "r'H" denotes the ray 
reflected off interface II which has traversed the medium 
and is incident from below on interface I. 

Rays that reflect off an interface have vector compo- 
nents in the x-z plane shifted in phase by 180°, which 
multiplies the x- and z-components of D and H by 
e 47r = — 1. Rays that traverse the uniaxial crystal develop 
a phase shift that depends on their angle, the thickness 
of the crystal, and the refractive index seen by that ray. 
The phase shifts differ for the ordinary and extraordi- 
nary rays, and are given by 



5' = h t cos 9' 
5" = h"t cos 9" ' 



(22) 



where h = n(l — ztan^) 1 / 2 is the complex refractive 
index, the imaginary part of which is often characterized 
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by the loss tangent, tan S, and encodes dielectric loss in 
the material. The real part is the refractive index given 
above. Rays that are transmitted through Interface I are 
incident on Interface II with a phase shift. This relates 
D t i and Dm as 



on Interface II to those on Interface I. Again, the 
relation can be written in matrix form as A// = 
(Df^Hl^D^-Hfj) = * n X with = P*j. The 
matrix P gives the effect on the fields of propagating 
through the material and is given by 



(23) 



D' UI = D' tI exp(*koS') 
D'iii = D&exp^") ' 

Likewise, the ray that reflects off interface II is phase 

shifted on its way to interface I, giving 



D 'rii = D' r , n eM-^oS') 

n" 

U rII — 



D';, n exp(-ik 5") 



(24) 



In the above, fco is the wave number for the wave in 
vacuum, equal to 2tt/Xq. For concision, we can break 
D and H into unknown total magnitudes multiplied by 
unit- vector directions that are known from Equations 20, 
21, and 12. The components of D and H transmitted 
from Interface I can be written out explicitly as 



/ exp(-A') 

exp(-A") 
exp(A') 

\ exp(A") j 

(28) 

where A' = «fc <5' and A" = ik 5" . The phases 5' and 8" 
are defined in Equation 22. 

We can find the relation we are seeking between A/ 
and Ajj by solving for X in the above equations and 
setting these equal to one another 



A/ = ^/(P^j)-^ 
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where use has been made of the relationship between 
the magnitudes of D and H, Equation 11. The complex 
components of the field unit vectors, such as D'^f 1 , are 
known and given in explicit form in Equation 30 below. 
All other vectors at the two interfaces, shown in Figure 
1, can be written in terms of the known transmitted field 
directions from interface I, as well as undetermined field 
amplitudes. The remaining field components at interface 
I can be written 



The matrix $/(P<f>/) _1 is the generalized transfer ma- 
trix for a uniaxial crystal with its optical axis at an angle 
X to its planar surface. It should be stressed that this ma- 
trix deals with total electric displacement and magnetic 
fields at the two boundaries of the crystal, which allows 
the matrix to take into account multiple reflections. Ex- 
plicit formulas for the field components transmitted from 
Interface I are 
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A system of four equations relating the total-field x- 
and y-components at Interface I and the individual ray 
components (tl and r'H) can be written in matrix form 
as A 7 = (Df , H V I ,D V I , -Hf) = $,A with 
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and X = (\D' tI \, |£>?j|, \D' r , H \, \D';, n \) 

A similar relation holds for Interface II, where Equa- 
tions 23 and 24 can be used to relate the fields 
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C. Matrix of an isotropic medium 

The above equations simplify considerably for the case 
of an isotropic dielectric. Specifically, s- and p-polarized 
waves no longer mix within the dielectric and the gener- 
alized transfer matrix becomes block diagonal. Taking x 
to be zero, 9i = 9 2 = 9, and n' = n" = n, Equations 30 
reduce to give D' tI = (0,1,0); D"j = (cos0, 0, - sin6»); 
H' tI = (-cos6»,0,sin6»); and H" T = (0,1,0). Carrying 
these field components through Equations 31, 28, and 
29 yields the generalized transfer matrix for an isotropic 
medium 
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The 2x2 matrices in the upper left and lower right of the 
isotropic generalized transfer matrix are those familiar 
from the literature for s- and p-polarized light, noting 
that the generalized transfer matrix as developed relates 
elements of D, whereas the treatments of isotropic media 
normally relate E = D/n 2 . 

4. Conclusion 

We have developed a matrix method capable of exact 
treatment of stratified optical systems involving bircfrin- 
gent crystals. This method allows straightforward calcu- 
lation of transmitted and reflected polarization states 
given an incident polarization state, for any frequency 
and angle of incidence. The generalized transfer matri- 
ces of individual layers can be multiplied together to give 
the matrix of a system composed of arbitrarily-many lay- 
ers. Because the generalized transfer matrix relates the 
total electric and magnetic fields at the two boundaries, 
multiple reflections between interfaces are automatically 
included in the treatment. 

Though the generalized transfer matrix method is able 
to treat general birefringent crystals, we have limited 
ourselves to calculating the important case of a uniax- 
ial crystal with its optic axis in the x-y plane. Full ex- 
pressions for this case have been given, as well as the 
reduction of this result to isotropic media. 
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